function Delta = AIA_step2(I,Phi)
    F = length(I);
    Delta = zeros(1,F);
    for f = 1:F
        R = AIA_step2_R(I{f,1},Phi);
        S = AIA_step2_S(Phi);
        X = S \ R;
        Delta(f) = atan2(-X(3),X(2));
    end
    Delta = Delta - Delta(1);
    Delta(Delta<0) = Delta(Delta<0) + 2*pi;
end

function R = AIA_step2_R(I,Phi)
    R11 = sum(I(:),"omitnan");
    R12 = I.*cos(Phi);R12 = sum(R12(:),"omitnan");
    R13 = I.*sin(Phi);R13 = sum(R13(:),"omitnan");
    R = [R11, R12, R13];
    R = R';
end

function S = AIA_step2_S(Phi)
%     [M, N] = size(Phi);S11 = M*N;
    S11 = ~isnan(Phi);S11 = sum(S11(:),"omitnan");
    S12 = cos(Phi);S12 = sum(S12(:),"omitnan");
    S13 = sin(Phi);S13 = sum(S13(:),"omitnan");
    S22 = cos(Phi).*cos(Phi);S22 = sum(S22(:),"omitnan");
    S23 = cos(Phi).*sin(Phi);S23 = sum(S23(:),"omitnan");
    S33 = sin(Phi).*sin(Phi);S33 = sum(S33(:),"omitnan");
    S = [S11, S12, S13;
         S12, S22, S23;
         S13, S23, S33];
end